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Abstract 

This  paper  addresses  the  development  of  homogenized  energy  models  which  characterize  the  fer- 
roelastic  switching  mechanisms  inherent  to  ferroelectric  materials  in  a  manner  suitable  for  subsequent 
transducer  and  control  design.  In  the  first  step  of  the  development,  we  construct  Helmholtz  and  Gibbs 
energy  relations  which  quantify  the  potential  and  electrostatic  energy  associated  with  90°  and  180° 
dipole  orientations.  Equilibrium  relations  appropriate  for  homogeneous  materials  in  the  absence  or 
presence  of  thermal  relaxation  are  respectively  determined  by  minimizing  the  Gibbs  energy  or  bal¬ 
ancing  the  Gibbs  and  relative  thermal  energies  using  Boltzmann  principles.  In  the  final  step  of  the 
development,  stochastic  homogenization  techniques  are  employed  to  construct  macroscopic  models 
suitable  for  nonhonrogeneous,  polycrystalline  compounds.  Attributes  and  limitations  of  the  charac¬ 
terization  framework  are  illustrated  through  comparison  with  experimental  PLZT  data. 
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1  Introduction 


Transducers  exploiting  ferroelectric  compounds  —  e.g.,  PZT,  PLZT  and  PMN  —  offer  unique  actuator 
and  sensor  capabilities  due  to  the  electromechanical  coupling  inherent  to  the  compounds.  A  fundamen¬ 
tal  attribute  of  these  materials,  when  operating  within  the  ferroelectric  regimes  typically  considered 
for  transducer  design,  is  the  presence  of  hysteresis  and  constitutive  nonlinearities  due  to  the  noncen- 
trosymmetric  structure  of  the  compounds.  These  nonlinearities  can  be  reduced  by  regulating  input 
fields  or  voltages,  employing  feedback  designs,  or  utilizing  charge  or  current-controlled  amplifiers,  in 
which  case,  linear  constitutive  relations  provide  reasonable  accuracy.  At  higher  drive  levels,  however, 
the  hysteresis  and  nonlinearities  associated  with  dipole  switching  must  be  accommodated  in  models, 
transducer  designs,  and  control  algorithms  to  achieve  design  specifications. 

As  detailed  in  Section  2,  hysteresis  and  constitutive  nonlinearities  in  ferroelectric  materials  can  be 
attributed  to  both  ferroelectric  and  ferroelastic  switching.  The  former  is  typically  due  to  field-induced 
90°  and  180°  dipole  switching  whereas  the  latter  is  associated  with  stress-induced  90°  switching.  Due 
to  the  electromechanical  coupling  inherent  to  compounds  such  as  PZT,  PLZT  and  PMN,  both  effects 
are  manifested  in  transducers  operating  in  high  drive  regimes. 

To  illustrate,  consider  the  PZT-based  THUNDER  (Thin  layered  UNimorph  Driver  and  sEnsoR) 
transducer  shown  in  Figure  1(a).  As  detailed  in  [24],  THUNDER  transducers  have  been  considered 
for  applications  ranging  from  high  speed  valve  design  to  shape  modification  in  an  airfoil.  Due  to  their 
construction,  they  are  capable  of  achieving  large  displacements  due  to  a  combination  of  robustness 
provided  by  the  metallic  backing  layer  and  curvature/stress  enhancement  achieved  during  fabrication. 
However,  the  prestresses  also  produce  the  asymmetric  hysteresis  behavior  shown  in  Figure  1(b)  which 
necessitates  the  use  of  models  which  incorporate  90°  ferroelastic  switching. 

Modeling  hierarchies  for  quantifying  ferroelastic  switching  can  be  characterized  by  the  level  at 
which  energy  or  phenomenological  principles  are  initiated,  since  essentially  all  the  frameworks  results 
in  macroscopic  models  or  constitutive  relations  commensurate  with  experimental  measurements  or 
transducer  design.  Micromechanical  models  are  typically  associated  with  energy  characterization  or 
thermodynamic  principles  developed  at  the  level  of  a  single  lattice  cell,  single  domain,  or  single  crystal. 
Subsequent  macroscopic  models  are  obtained  through  brute  force  computations  or  various  homoge¬ 
nization  techniques.  Material  characterization  via  large  scale  computations  can  be  used  to  quantify 
fundamental  mechanisms  whereas  the  low-order  homogenized  models  facilitate  implementation.  Ad¬ 
ditionally,  the  effective  parameters  resulting  from  homogenization  can  phenomenologically  quantify 
unknown  physics  or  model  uncertainty,  thus  yielding  models  that  are  accurate  for  a  range  of  operating 
conditions.  Additional  details  regarding  micromechanical  models  for  materials  can  be  found  in  the 
summaries  [10,  12,  24]  and  references  [11,  16,  17]. 


Figure  1:  (a)  THUNDER  transducer,  and  (b)  stress-dependent  electromechanical  behavior. 
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A  second  technique  for  constructing  macroscopic  models  is  to  employ  thermodynamic  or  phe¬ 
nomenological  principles  directly  at  the  macroscopic  level  to  characterize  measured  phenomena.  In 
the  former  category,  states  such  as  remanent  polarization  or  remanent  strain  re-identified  as  internal 
variables  and  the  first  and  second  laws  of  thermodynamics  are  invoked  to  characterize  material  be¬ 
havior  [8,  33]  —  the  reader  is  referred  to  Landis  [10]  for  a  summary  of  constitutive  theories  resulting 
from  reversible  and  irreversible  thermodynamics.  Purely  phenomenological  approaches  range  from 
the  rheology  theory  of  [32]  to  Preisach  models  [6,  21].  We  note  that  Preisach  theory  for  ferroelectric 
materials  is  significantly  less  mature  than  its  ferromagnetic  counterpart  and  presently  is  focused  on 
ferroelectric  rather  than  ferroelastic  switching  mechanisms. 

The  present  framework  combines  energy  principles  at  the  lattice  level,  theory  of  thermally  activated 
processes,  and  stochastic  homogenization  techniques  to  characterize  hysteresis  due  to  ferroelectric  and 
ferroelastic  switching  in  a  manner  which  facilitates  material  characterization,  transducer  design,  and 
model-based  control  design.  In  the  first  step  of  the  development,  Helmholtz  and  Gibbs  energy  relations 
are  constructed  at  the  lattice  level  to  quantify  the  internal  and  electrostatic  energy  associated  with  90° 
and  180°  dipole  orientations.  To  characterize  regimes  in  which  thermal  activation  is  significant,  the 
Gibbs  and  relative  thermal  energies  are  balanced  through  Boltzmann  theory  to  provide  equilibrium 
relations  quantifying  local  strains  and  polarizations  as  a  function  of  input  stresses  and  fields.  As  shown 
in  [24,  31] ,  these  local  relations  reduce  to  minima  of  the  Gibbs  energy  in  the  limit  of  negligible  thermal 
activation  —  enforcement  of  these  minimization  criteria  significantly  improve  the  efficiency  of  the 
model  and  inversion  process  employed  for  linear  control  design  in  regimes  where  relaxation  processes 
are  negligible.  This  provides  a  model  for  homogeneous,  single  crystal  compounds  with  uniform  effective 
fields.  To  incorporate  the  effects  of  poly  crystallinity,  material  nonhomogeneities,  and  variable  effective 
fields,  parameters  such  as  the  local  coercive  fields  are  considered  to  be  manifestations  of  underlying 
distributions  rather  than  constants.  Homogenization  in  this  manner  yields  low-order  models  posed  in 
terms  of  effective  parameters  that  are  efficient  to  implement  in  optimization  or  control  algorithms. 

An  important  attribute  of  the  framework  is  the  flexibility  it  provides  for  characterizing  hystere¬ 
sis  and  constitutive  nonlinearities  in  a  broad  range  of  ferroelectric,  ferromagnetic  and  ferroelastic 
materials.  This  combination  of  energy  analysis  and  theory  of  thermally  activated  processes  had  its 
genesis  in  the  SMA  models  developed  by  Muller,  Achenbach  and  Seelecke  [1,  23].  The  incorporation 
of  stochastic  homogenization  techniques  and  extension  to  ferroelectric  [27,  31],  ferromagnetic  [25,  26], 
and  polycrystalline  SMA  compounds  [14,  15,  22]  has  led  to  the  formulation  of  a  general  framework 
for  characterizing  hysteresis  in  ferroic  compounds  [24,  29,  30].  The  theory  presented  here  extends 
this  framework  and  the  work  presented  in  [9]  by  incorporating  stress  and  field-induced  90°  switch¬ 
ing  as  motivated  by  both  fundamental  material  considerations  and  design  criteria  associated  with 
stress-dependence  as  illustrated  in  Figure  1(b). 

Physical  mechanisms  associated  with  ferroelectric  and  ferroelastic  switching  are  summarized  in  Sec¬ 
tion  2  to  motivate  90°  and  180°  switching  mechanisms  that  need  to  be  incorporated  in  energy  relations. 
In  Section  3,  we  summarize  the  construction  of  a  Gibbs  energy  functional  based  on  Landau-Devonshire 
principles.  Minimization  of  this  functional  quantifies  homogeneous,  single  crystal  behavior  in  the  ab¬ 
sence  of  relaxation  processes.  Whereas  efficient  to  formulate,  the  Landau-Devonshire  functional  is 
overly  restrictive  and  often  unstable  to  implement  due  to  its  reliance  on  high-order  polynomials.  To 
alleviate  these  restrictions,  2-D  piecewise  quadratic  energy  functionals,  motivated  by  the  Landau- 
Devonshire  formulation,  are  employed  in  Section  4  to  provide  local  models  characterizing  ferroelectric 
and  ferroelastic  switching.  In  Section  5,  stochastic  homogenization  techniques  are  employed  to  con¬ 
struct  a  macroscopic  model  for  polycrystalline,  nonhomogeneous  materials  and  in  Section  6,  properties 
of  the  model  are  illustrated  through  comparison  and  prediction  of  experimental  PLZT  data. 
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2  Ferroelectric  and  Ferroelastic  Switching  Mechanisms 


To  motivate  mechanisms  which  must  be  incorporated  in  energy  formulations,  we  summarize  ferro¬ 
electric  and  ferroelastic  switching  mechanisms  in  ferroelectric  compounds.  We  will  focus  primarily 
on  Pb(Zr,Ti)03  (lead  zirconate-titanate  or  PZT),  which  is  comprised  of  PbZr^-iOs  (lead  zirconate) 
and  PbTi_r03  (lead  titanate)  where  x  is  chosen  to  optimize  electromechanical  coupling,  and  lanthanum 
doped  PZT  (PLZT).  For  temperatures  above  the  Curie  point  Tc,  the  structures  of  PbTiC>3  and  PbZrC>3 
are  cubic  whereas  for  T  <  Tc,  the  structure  of  PbTiC>3  is  tetragonal  and  PbZrC>3  is  orthorhombic  [24]. 
We  illustrate  the  switching  mechanisms  in  the  context  of  the  paraelectric  cubic  and  the  ferroelectric 
tetragonal  structure  of  lead  titanate  in  Figure  2  and  note  the  analogous  behavior  is  observed  for  the 
orthorhombic  structure  of  lead  zirconate.  To  simplify  the  model  while  retaining  the  structure  required 
for  stress-induced  switching,  we  consider  a  2-D  polarization  P  =  (Pi,  P3)  having  the  orientation  shown 
in  Figure  2(a). 

Ferroelectric  switching  is  induced  by  the  application  of  an  electric  field  E  that  is  larger  in  magnitude 
than  the  coercive  field  Ec.  For  PbTiC>3,  this  causes  the  central  Ti+4  ion  to  relocate  to  a  new  equilibrium 
position,  resulting  in  a  180°  change  in  polarization  that  is  parallel  to  the  applied  field  as  depicted  in 
Figure  3(a).  Ferroelastic  switching  is  caused  by  the  application  of  a  stress  cr  that  is  larger  in  magnitude 
than  the  coercive  stress  ac  producing  a  90°  change  in  polarization  that  is  perpendicular  to  the  applied 
stress  as  illustrated  in  Figure  3(b).  The  ferroelectric  and  ferroelastic  switching  mechanisms  cause 
a  hysteretic  relationship  between  input  fields  E  and  cr  and  output  polarization  P  and  strains  e. 
Experimental  data  illustrating  typical  ferroelectric  and  ferroelastic  response  of  the  soft  piezoelectric 
ceramic  PLZT,  can  be  found  in  [13]. 

Figures  4(a)  and  4(b)  illustrate  the  relationship  between  an  externally  applied  field  and  the  po¬ 
larization  and  strains  induced  in  a  soft  PLZT  compound.  At  point  A,  the  electric  field  is  sufficiently 
strong  so  that  all  the  dipoles  form  one  domain  that  is  aligned  in  the  direction  of  the  applied  field.  As 
the  field  is  decreased,  it  approaches  the  negative  coercive  field  in  the  region  around  point  B  where 
180°  switching  commences.  Additionally,  this  often  includes  90°  switching  as  indicated  by  the  pres¬ 
ence  of  a  negative  strain  at  point  B  as  depicted  in  Figure  4(b).  At  point  C,  all  the  dipoles  have 
switched  and  again  form  one  domain  that  is  aligned  in  the  direction  of  the  electric  field.  At  point  C, 
the  polarization  is  opposite  to  that  at  point  A  whereas  the  strains  have  the  same  value.  As  the  field 
is  increased,  90°  switching  occurs  at  point  D  and  rapidly  continues  back  to  point  A  where  the  full 
180°  switch  has  occurred  and  the  dipoles  are  again  aligned  with  the  applied  electric  field. 


Figure  2:  (a)  Coordinate  system,  (b)  high  temperature  paraelectric  cubic  form  of  PbTiC>3,  and  (c)  low 
temperature  ferroelectric  tetragonal  form  of  lead  titanate  and  resulting  spontaneous  polarization  Pq. 
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(a) 


(b) 


Figure  3:  (a)  Ferroelectric  180°  switch  in  spontaneous  polarization  Po  induced  by  an  applied  electric 
field,  and  (b)  ferroelastic  90°  switch  induced  by  an  applied  stress. 


Figures  4(c)  and  4(d)  illustrate  the  relationship  between  an  externally  applied  stress  and  the 
polarization  and  strains  induced  in  the  soft  PLZT.  At  point  F,  the  dipoles  are  aligned  in  the  positive 
3-direction  and  the  material  acts  as  one  domain.  As  the  compressive  stress  is  increased  in  magnitude, 
it  approaches  the  coercive  stress.  In  the  region  around  point  G,  90°  switching  occurs  and  the  dipoles 
begin  to  align  perpendicular  to  the  direction  of  the  applied  stress.  This  is  indicated  by  the  presence 
of  a  negative  strain  at  point  G  in  Figure  4(d).  As  the  stress  is  reduced  in  magnitude,  the  material 
remains  poled  perpendicular  to  the  applied  stress  resulting  in  a  decrease  in  the  polarization  in  the 
3-direction. 

3  Single  Crystal  Model:  Landau— Devonshire  Energy  Relation 

To  motivate  issues  which  must  be  addressed  when  constructing  energy  relations  and  obtaining  condi¬ 
tions  necessary  for  equilibrium  states,  in  the  absence  of  thermal  excitation,  we  first  employ  Landau- 
Devonshire  principles  to  construct  a  local  model.  This  approach  has  the  advantage  of  simple  energy 
representations  but  is  overly  restrictive  due  to  the  requisite  high-order  polynomials.  These  restric¬ 
tions  are  addressed  in  Section  4  where  the  Landau  Devonshire  theory  is  used  to  motivate  piecewise 
energy  definitions.  Using  the  theory  of  Landau  and  Devonshire,  the  Helmholtz  energy  is  expressed  as 
a  truncated  power  series  in  term  of  the  polarization  with  coefficients  chosen  to  ensure  measured  ma¬ 
terial  properties.  As  detailed  in  [24],  the  behavior  of  materials  exhibiting  first  and  second-order  phase 
transitions  can  be  quantified  by  the  retention  of  Ath  and  6th  order  polarization  terms  of  the  Helmholtz 
energy.  Lead  titanate  and  lead  zirconate  both  exhibit  first-order  transitions  near  the  transition  tem¬ 
perature.  However,  we  shall  consider  operating  regimes  T  «  Tc,  where  the  materials  are  actually 
employed  for  transduction;  which  necessitates  retention  of  polarization  terms  only  up  to  4<ft-  order. 


Figure  4:  (a)  Hysteretic  field-polarization  relation  for  bulk  PLZT,  (b)  field-strain  behavior  of  PLZT, 
(c)  stress-polarization  relation  for  PLZT,  and  (d)  stress-strain  behavior  of  PLZT. 
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To  model  the  internal  free-energy  of  a  single  lead  titanate  crystal,  we  consider  electric  field  E  = 
(Ei,  £3)  and  stress  <j  =  ( <71,(73)  inputs  and  polarization  P  =  (Pi ,  P3)  and  strain  e  =  (£1,63)  outputs. 
The  polarization  component  of  the  Helmholtz  free-energy  i[’p  is  taken  to  be  the  Ath  order  polynomial 

V>p(P)  =  oqPj2  +  a3Pz  +  au  Pf  +  «33  Pi  +  a3iP32P12. 

The  coefficients  are  chosen  so  that  at]  >  0  and  a*  <  0  below  the  Curie  point  and  can  be  related  to 
physical  properties  of  ferroelectric  compounds  such  as  the  renranence  polarization  Pp  and  coercive 
field  Ec  —  e.g.,  see  (8)  and  (9).  The  electromechanical  coupling  component  is  given  by 

As{ P,  e)  =  -aieiPi  -  a3£3P3  -  ai3£iP3  -  a3i£3Pi  -  q\£\P‘f  -  g3£3P32  -  <?i3£iP32  -  <73i£3P2  (1) 

where  are  piezoelectric  coupling  coefficients  and  qij  are  electrostrictive  coupling  coefficients.  When 
shear  effects  are  neglected,  the  elastic  free-energy  is 

^elis)  =  2^l£l  +  2^3£^  (2) 

where  Yz  are  the  components  of  the  Young’s  modulus.  The  total  Helmholtz  free-energy  is  then  given 
by 

■0(P,  e)  =  V+K P)  +  1pes( P,  e)  +  ^eK£)- 

The  work  due  to  an  externally  applied  electric  field  and  an  applied  stress  is  incorporated  by  employing 
a  Gibbs  energy  of  the  form 


G(E,  P,  cr,  e)  =  ^(P,  e)  —  E  •  P  —  <T  ■  E. 


(3) 


As  detailed  in  [3,  24],  necessary  conditions  for  minimizing  the  Gibbs  energy  in  the  absence  of 
thermal  relaxation  mechanisms  are 


dG_  dG_  _  dG_ 
dei  ’  de3  ’  dPi 


=  0, 


dG 

dP3 


=  0. 


These  conditions  specify  how  the  polarization  and  strains  change  to  minimize  the  internal  energy  when 
an  external  force  is  applied  and  thermal  activation  is  negligible.  The  condition  =  0  implies  that 


£1  —  Y\  1  (<7i  +  ai-Pi  +  ai3P3  +  q\ Pf  +  qi3P3)  (4) 

and  =  0  implies  that 

£3  =  Y3  1  (<73  +  a3P3  +  a3\P\  +  q3P$  +  q3iPi)  .  (5) 


These  relations  allow  us  to  determine  explicit  expressions  for  the  strains.  In  a  manner  analogous  to 
that  employed  in  [19]  for  a  stress-free  state,  we  substitute  equations  (4)  and  (5)  into  the  Gibbs  energy 
(3)  to  directly  couple  the  stress  and  the  polarization.  This  yields 

G(E,  P,  <r)  =  ^(P,tr)  -E-P 

where  cr)  is  the  fourth  order  polynomial 

£(P,  <?)  =  7i  Pf  +  72  Pi  +  73  P$P?  +  74^i3  +  75-F3  +  ^p?Pt  +  yrPiPi 

+78  P1  +  79-ff  +  710-P3-P1  +  711 -Pi  +  712P3  +  713- 
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(a) 


Figure  5:  Gibbs  energy  landscape,  (a)  Coordinate  system  of  energy  landscape,  and  (b)  contour  plot 
of  Gibbs  energy  landscape. 


The  Gibbs  energy  at  zero  stress  and  zero  electric  field  has  four  local  minima  each  corresponding  to 
four  energetically  favorable  polarization  states  EP\  and  ±P->  and  is  illustrated  in  Figure  5.  Given  an 
input  pair  (E,  <r),  the  local  polarization  in  the  absence  of  thermal  excitation  can  then  be  calculated 
by  locally  minimizing  the  Gibbs  energy  to  obtain 

(P)  =  argnrin  G(E,  P,  cr).  (6) 

The  ferroelastic  model  thus  quantifies  the  (Pi,  P3)  pair  that  minimizes  the  total  energy  of  the  system 
and  the  corresponding  strain  pair  (£1,63)  is  designated  by  equations  (4)  and  (5). 

As  noted  previously,  several  of  the  coefficients  in  the  ferroelastic  model  are  related  to  physical 
material  properties,  thus  allowing  a  means  of  determining  their  values.  The  remanent  polarization  Pr 
and  coercive  field  Ec  are  computed  by  considering  a  stress-free  Gibbs  energy 


G(E,P)  =^(P)  -E-P 

and  noting  that  the  conditions  Pi  =  0  and  =  0  imply 

4a33-P|  +  2«3P3  —  £3  =  0.  (7) 


As  outlined  in  [34],  equation  (7)  has  three  real  roots,  of  which  at  least  two  are  equal  if  the  discriminant 
is  0.  The  coercive  field  value  for  which  this  property  is  satisfied  is 


Ec  = 


(8) 


When  £3  =  0,  the  remanent  polarization  determined  from  solving  (7)  is  given  by 


pr = m  <9) 

Equations  (8)  and  (9)  provide  a  system  of  two  equations  and  two  unknowns  sufficient  to  solve  for  a 3 
and  «33  given  values  for  Ec  and  Pr. 

A  similar  process  may  be  applied  to  find  aq  and  an  if  the  energy  landscape  is  known  to  be 
nonsymmetric.  Whereas  some  of  the  coefficents  of  this  Landau  Devonshire  formulation  are  related  to 
physical  quantities,  it  is  difficult  to  relate  all  parameters  with  material  properties.  Furthermore,  the 
model  proves  overly  restrictive  to  implement  in  various  regimes  due  to  the  high  order  polymonials. 
In  the  next  sections,  we  present  two  piecewise-defined  ferroelastic  switching  models  that  accurately 
describe  the  physics  and  yield  methods  for  predicting  material  parameters. 
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4  Single  Crystal  Model:  Piecewise  Quadratic  Energy  Relations 


The  high-order  polynomials  in  the  Landau-Devonshire  theory  of  Section  3  can  induce  unstable  be¬ 
havior  through  small  modifications  of  the  coefficients.  To  facilitate  model  implementation,  improve 
computational  efficiency  and  increase  the  flexibility,  we  employ  piecewise  quadratic  energy  function¬ 
als  which  are  approximations  of  the  energy  functionals  based  on  Landau-Devonshire  principles.  The 
model  coefficients  directly  relate  to  measured  material  properties  allowing  a  means  for  parameter  iden¬ 
tification  and  estimation.  The  low-order  piecewise  polynomial  models  also  facilitate  transducer  and 
control  design. 

4.1  2-D  Ferroelectric  Switching  Model 

We  summarize  first  a  2-D  ferroelectric  switching  model  which  characterizes  the  hysteretic  field- 
polarization  and  field-strain  relations  depicted  in  Figures  4(a)  and  4(b)  and  allows  90°  polarization 
switching  to  occur  when  energetically  favorable.  However,  the  ferroelectric  model  does  not  characterize 
the  stress-induced  switching  mechanisms  resulting  in  the  stress-polarization  and  stress-strain  relations 
depicted  in  Figures  4(c)  and  4(d)  and  should  be  limited  to  low  stress  regimes.  The  extension  of  the 
model  to  incorporate  ferroelastic  switching  is  addressed  in  Section  4.2. 


4.1.1  Helmholtz  and  Gibbs  Energy  Relations 

For  fixed  temperatures,  we  expand  on  the  theory  of  [31]  where  it  is  illustrated  that  a  reasonable 
expression  for  the  Helmholtz  energy  is  a  piecewise  quadratic  relation.  In  the  context  of  the  model 
presented  in  Section  3,  this  is  obtained  by  retaining  2nd  order  terms  of  the  Taylor  approximations  for 
the  4th  order  polynomial  expressions  about  the  equilibria.  To  account  for  the  180°  and  90°  polarization 
switching,  we  allow  four  possible  polarization  orientations  (P^,  P±)  as  depicted  in  Figure  5.  We  define 
piecewise  quadratic  polynomials  along  the  P\  and  P3  directions  to  obtain  the  functionals 


Vh(Pi) 


MP3) 


\rti  (Pi  +  Pi-R)2 

,  Pi  <  -Pi/ 

\ili  (Pi  -  Pir)2 

,  Pi  >  Pi/ 

\fli  ( P\i  ~  Pir )  | 

[ht  -  pi«) 

,  IP1I<P1/ 

5%  (P3  +  P3R.)2 

)  P3  <  -P3I 

5%  (P3  -  P3R.)2 

,  P3  >  P3I 

2?73  (P31  ~  P3R)  ( 

%  - p™) 

,  IP}]  <  P3I 

where  P\j  and  P3/  denote  the  positive  inflection  points  in  the  P\  and  P3  directions  and  P\r  and  P3# 
denote  the  polarization  at  which  the  positive  minimum  occurs  in  the  Pi  and  P3  directions.  The  2-D 
Helmholtz  energy  is  then  defined  by 


Mp)  =  Vh(Pi)  +  MP3)  +  cPlPl 


(10) 


The  resulting  Gibbs  free  energy,  as  illustrated  in  Figure  5,  quantifying  the  change  in  the  energy 
landscape  due  to  an  applied  field,  is  provided  by  the  inclusion  of  the  electrostatic  work  term  and 
results  in 


G(E,P)=Y>p(P)-E-P. 


(11) 
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4.1.2  Switching  in  the  Absence  of  Thermal  Activation 


In  the  absence  of  thermal  activation,  ferroelectric  dipole  switching  occurs  when  the  magnitude  of 
the  applied  electric  field  is  greater  than  that  of  the  coercive  field  Ec.  From  an  energy  perspective, 
the  polarization  switching  occurs  when  the  applied  electric  field  is  sufficiently  large  to  eliminate  local 
minima  in  the  Gibbs  energy.  The  local  polarization  (P)  of  the  single  crystal  can  be  calculated  by 
solving  the  necessary  conditions 


dG 


=  0, 


dG 


=  0. 


dPi  ’  0P3 

Alternatively,  given  an  input  field  E,  the  local  polarization  can  be  found  by  locally  minimizing  the 
Gibbs  energy, 


(P)  =  argnhn  G(E,  P) 


(12) 


where  G  is  defined  by  equation  (11).  Since  the  single  crystal  has  a  dipole  structure  that  will  change 
orientation  to  minimize  the  internal  energy,  the  order  in  which  the  minima  are  eliminated  determines 
the  type  of  ferroelectric  switching  that  occurs. 


Ferroelectric  90°  Switching 

A  ferroelectric  90°  switch  is  illustrated  in  Figure  6(a)  where  an  applied  field  has  caused  the  P.^  min¬ 
ima  to  vanish  while  the  P±  minima  remain.  If  the  single  crystal  was  oriented  in  the  P 3”  direction, 
this  would  result  in  a  90°  switch  to  one  of  the  neighboring  wells  in  the  Pi  direction.  Given  an  ap¬ 
plied  electric  field  E3,  we  can  calculate  the  field  necessary  to  induce  switching  by  first  solving  the 
necessary  condition  =  0.  We  are  calculating  the  value  of  the  electric  field  that  will  make  the 
Pg~  minima  disappear,  or  equivalently  the  electric  field  that  will  make  P3  =  P3/.  This  occurs  when 
Ec  =  V3  (PiR  —  P31 )  +  cPiP3j.  However,  the  second  necessary  condition  =  0  requires  that  Pi  =  0 
since  we  assume  the  applied  field  has  no  E\  component.  This  results  in 

Ec  =  V3  {P3R  ~  P31)  •  (13) 

The  electric  field  required  to  eliminate  the  minima  in  the  90°  orientation  is  denoted  by  P®°.  For  a 
ferroelectric  90°  switch  to  occur,  >  Ec. 


Figure  6:  (a)  Elimination  of  180°  minimum  through  application  of  an  applied  field  greater  than  Ec, 
and  (b)  90°  minima  can  disappear  before  the  180°  minima  disappears  if  <  Ec. 


Ferroelectric  180°  Switching 

As  detailed  in  Section  2,  180°  polarization  switching  due  to  an  applied  electric  field  can  result  from 
a  direct  180°  switch  (e.g.,  P3“  — >  P 3")  or  a  pair  of  90°  switches  (e.g.,  P%  — >  P±  — >  P%)-  The  180° 
switching  mechanism  which  consists  of  a  pair  of  90°  switches  occurs  when  |  £^3 1  >  E®°  >  Ec.  The 
direct  180°  switch  is  illustrated  in  Figure  6(b)  where  Ec  >  E®°  and  the  applied  field  has  caused  the 
P^  and  P±  minima  to  vanish.  If  the  single  crystal  was  oriented  in  the  P^  direction,  this  would  result 
in  a  direct  180°  switch  to  the  well  in  the  P 3"  direction. 

4.1.3  Switching  in  the  Presence  of  Thermal  Activation 

In  the  model  developed  in  Section  4.1.2,  dipole  switching  occurred  only  when  the  magnitude  of  the 
applied  electric  field  exceeded  that  of  the  coercive  field  value  of  the  single  crystal.  However,  mechanisms 
such  as  excitation  from  thermal  effects  can  induce  switching  before  local  minima  are  eliminated.  To 
model  this  phenomena,  we  incorporate  the  effects  of  thermal  activation.  We  continue  to  assume  that 
the  material  is  homogeneous  and  that  each  dipole  in  this  collection  has  the  same  energy  landscape 
and  are  subject  to  the  same  switching  mechanisms.  This  allows  us  to  model  the  switching  that  occurs 
by  evolving  the  fraction  of  dipoles  in  each  allowed  orientation  by  a  population  model  that  is  similar 
to  what  is  derived  in  [31]. 

As  detailed  in  [24],  the  Gibbs  energy  and  the  relative  thermal  energy  ^  is  balanced  through  the 
Boltzmann  relation 

p(G(  P))  =  Ce~G(P)V/kT.  (14) 

Here  V  denotes  a  representative  control  volume,  k  is  Boltzmann’s  constant  and  T  is  the  temperature. 
The  fraction  of  dipoles  in  each  allowed  dipole  orientation  is  given  by  x$,  x$,  xf,  and  xj~.  The  average 
polarization  associated  with  each  allowed  dipole  orientation  (the  minima  of  the  Gibbs  energy)  is 
denoted  by  (P3  ),  (Pg  ),  (P^),  and  (P^).  It  should  be  noted  that  these  are  vector  quantities  and 
have  polarization  values  in  both  the  Pi  and  P3  directions.  Conservation  of  the  total  number  of  dipoles 
yields 


x3  +  x3  +  4  +  X1  =  1. 

(15) 

For  a  homogeneous  material,  the  components  of  the  local  polarization 

(P)  =  (Pl,P3) 

(16) 

are  given  by 

^3  =  4  <P3+  >3  +  X3  < P3 ~ >3  +  xf  (P+  >3  +  aq  (Pf  >3 

Pi  =  4  (P4i  +  *3  <^3-)!  +  4  (pi)i  +  (pi)i 

where  the  vector  components  of  the  minima  are  denoted  by  the  subscripts  1  or  3. 

The  likelihood  of  switching  from  the  Pj-  orientation  (south)  to  the  P^  (east)  is  specified  by  pse 
and  is  calculated  by 

=  7 _ 1 _ 

PS£  T  e— 7(G(P^in)— G(P-rrier))  _  j 

where  the  relaxation  time  r  is  the  reciprocal  of  the  frequency  at  which  dipoles  attempt  a  switch  and 
7  =  -up-  Here,  Plearrier  denotes  the  location  of  the  saddle  point  from  the  Pf"  orientation  (south)  to 
P3"  (east)  and  G  (Psbearrier)  is  the  Gibbs  energy  evaluated  at  that  point.  Similarly,  Psrnm  denotes  the 
location  of  the  south  minima  and  G  (P(rem)  is  the  Gibbs  energy  evaluated  at  that  point. 
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The  likelihood  of  switching  from  the  P ^  orientation  {north)  to  the  P^  {east)  is  specified  by  pne 
and  is  calculated  by 

y  =  7 _ 1 _ 

Pn£  T  e— 7(G(P"  in)~G(PZr rier))  _  l 

where  P barrier  denotes  the  location  of  the  saddle  point  from  the  P+  orientation  {north)  to  P%  {east) 
and  G  {Pfarrier)  is  the  Gibbs  energy  evaluated  at  that  point.  Similarly,  P^in  denotes  the  location  of 
the  north  minima  and  G  (P”„n)  is  the  Gibbs  energy  evaluated  at  that  point.  The  remaining  likelihoods 
follow  similarly  and  have  the  notation  summarized  in  Table  1. 

Finally,  the  evolution  of  the  dipole  fractions  due  to  thermal  activation  is  governed  by  the  system 
of  differential  equations 

ry  I  -  ry\  ry  I  ry\  ry  I  ry\  ry  I  ry\  ry  I 

^3  —  yse^ i  i  yne^i  yes^s  t)endj2> 

% 3  =  Psw%i  H-  Pnw% i  Pws% 3  Pwn% 3 

ry  I  -  ry\  ry  I  ry\  ry »  I  _  ry\  ry *  I  _  ry~\  ry  I 

—  ywn^ 3  “t  yen*1' 3  ynw^i  ynedj  1 

■7 1  —  Pws% 3  T  Pes'K 3  Pse%\  • 

The  dimension  of  this  system  of  first  order  ordinary  differential  equations  can  be  reduced  by  using 
(15)  and  solved  quickly  and  with  sufficient  accuracy  using  an  implicit  method  such  as  backward  Euler 
method.  To  compute  the  transition  likelihoods,  values  of  the  minima  and  the  saddle  points  of  the 
energy  landscape  are  computed  using  Newton’s  method  with  line  search  to  solve  for  the  location  where 
the  gradient  of  the  Gibbs  energy  is  zero. 

4.2  2-D  Ferroelastic  Switching  Model 

Here  we  extend  the  framework  of  Section  4.1  to  construct  a  ferroelastic  switching  model  which  quanti¬ 
fies  the  hysteretic  field-polarization  and  field-strain  relationships  depicted  in  Figures  4(a)  and  4(b)  and 
also  characterizes  the  hysteretic  stress-polarization  and  stress-strain  relationships  depicted  in  Figures 
4(c)  and  4(d).  The  ferroelastic  model  extends  the  phenomenological  Landau-Devonshire  energy  model 
developed  in  [2]  and  summarized  in  Section  3  through  the  construction  of  piecewise  quadratic  energy 
functionals. 

4.2.1  Helmholtz  and  Gibbs  Energy  Relations 

To  model  stress-induced  switching  in  ferroelectric  materials,  we  again  consider  electric  field  E  = 
{E\,Ei)  and  stress  a  =  (<71,03)  inputs  and  polarization  P  =  ( Pi,P 3)  and  strain  e  =  (£1,63)  outputs. 


State  Transition 

notation 

Pi  ~ 

-4-  P+ 

-  <3 

Pse 

Pi  - 

-Pi 

Psw 

Pt ~ 

-Pt 

Pne 

Pt- 

-Pi 

Pnw 

Pi  ~ 

-Pt 

Pwn 

Pi~ 

-Pi 

Pws 

p+  _ 

-Pt 

Pen 

p+  _ 

-Pi 

Pws 

Table  1:  2-D  transition  likelihoods. 
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The  polarization  component  of  the  Helmholtz  free-energy  ?/>p  is  defined  piecewise  by  (10)  and  is 
a  function  of  the  four  possible  polarization  orientations  (P±,P±).  The  electromechanical  coupling 
component  ipes  is  given  by  (1).  Ignoring  shear  effects,  the  elastic  free-energy  ipei  is  given  by  (2).  The 
total  Helmholtz  free-energy  is  then 


^(P,  e)  =  ^p(P)  +  1pes(  P,  e)  +  'Ipelis) 

and  the  Gibbs  energy  is 

G(E,  P,  a,  e)  =  ^(P,  e)  —  E  •  P  —  cr  •  e.  (17) 

As  detailed  in  Section  3  in  the  context  of  the  Landau  Devonshire  functional,  the  necessary  con¬ 
ditions  allow  us  to  determine  explicit  expressions  for  the  strains.  This  is  a  direct  result  of  neglecting 
the  shear  terms  in  i)>ei  which  allows  the  strains  to  decouple.  Substituting  equations  for  the  strains  into 
the  Gibbs  energy  (17)  yields 

G(E,P,<x)  =  ^(P,<r)  — E-P  (18) 

where  i/>(P,  cr)  is  a  fourth-order  piecewise  polynomial  defined  by 

£(P,  Or)  =  71  Pi4  +  72P34  +  73  PfPl  +  74P?  +  75P|  +  76^3  +  77AP32 


+78  Pi  +  79P32  +  710P3P1  +  711P1  +  712P3  +  713, 
with  the  7 j  coefficients  defined  in  Appendix  A. 


4.2.2  Local  Polarization  Model 


In  the  absence  of  thermal  effects,  the  local  polarization  (P)  with  an  applied  stress  can  be  determined 
by  solving  the  necessary  conditions 


dG 

dP\ 


=  0, 


dG 

dP3 


=  0 


or  by  minimizing  the  Gibbs  energy 


(P)  =  argnun  G(E,  P,  cr) 


(19) 


where  the  Gibbs  energy  is  defined  by  equation  (18).  To  incorporate  thermal  relaxation,  the  local 
average  polarization  can  be  calculated  by  using  the  thermal  evolution  model  developed  in  Section  4.1.3. 

To  calculate  the  value  of  the  coercive  field  Ec  for  the  stress-dependent  ferroelastic  model,  we  apply 
the  necessary  condition  J^r  =  0  and  solve  for  the  value  of  the  electric  field  that  eliminates  the  minima 
in  the  Gibbs  energy.  This  yields 


Ec{o -3)  = 


2(73p37  -  «3 


*3 

Pjaj  +  *3013 


(03)  + 


2Y3qf3  +  2Yxq\ 


Y>,Yi 

P3I  +  V3{P3R  ~  P3I ) • 


p3 

■‘■SI 


3y3ai3gi3  +  3Yia3g3 
T3P1 


P31  + 


(20) 


It  is  clear  from  equation  (20)  that  the  coercive  field  is  stress-dependent.  A  ferroelastic  switch  will 
occur  when  the  applied  stress  is  greater  in  magnitude  than  the  coercive  stress  ac.  To  calculate  ac,  we 
set  Ec  =  0  in  equation  (20)  resulting  in 


(Jr  = 


Y3 


«3  —  2(73+3/ 


2*3  gi3  +  2Lj(732 

T3P1 


p3 

3/ 


3*3Ql3gl3  +  31W3 

*3*1 


Pi,  + 


Y\  di 


Y3a2 
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T3** 


P37  +  V3{P3R  -  P31) 


(21) 
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4.3  1-D  Ferroelastic  Switching  Model 

To  provide  a  simplified  model  which  facilitates  real-time  implementation  for  transducer  and  control 
design,  we  present  here  a  1-D  model  that  incorporates  both  90°  ferroelectric  and  ferroelastic  polar¬ 
ization  switching  mechanisms.  For  the  1-D  input  fields  considered  in  the  examples  of  Section  6,  the 
1-D  ferroelastic  model  has  a  90°  switching  behavior  similar  to  the  2-D  ferroelastic  model  and  is  more 
computationally  efficient.  We  note  that  aspects  of  the  functionals  are  similar  to  those  employed  for 
SMA  undergoing  austinite  -  martinsite  phase  transformations  and  the  reader  is  referred  to  [14,  24,  29] 
for  details  illustrating  properties  of  the  SMA  relations. 


4.3.1  Helmholtz  and  Gibbs  Energy  Relations 


We  consider  electric  field  and  stress  inputs  ( E,a )  and  polarization  and  strain  outputs  (Re).  The 
polarization  has  three  allowed  dipole  states  P_,  P+  and  Tgo  as  shown  in  Figure  7.  The  values  P_  and 
P+  correspond  to  the  ±P3  orientation  in  the  2-D  model  while  Pgo  represents  both  ±Pi .  We  define  the 
polarization  component  of  the  Helmholtz  free  energy  to  be, 


MP)  =  { 


where 


P  — 

m  — 


V  {Pi  ~  Pr )  -  V2P901 


1{P  +  Pr? 

Cg 

1 

VI 

Cg 

f  ( P  +  Pm)2  +  P 

,  —Pi  <  P  <  -P901 

f  (P)2  +  A 

,  \P  <  P901 

f  (P-Pm)2  +  (3 

,  P901  <  P  <  Pi 

1{p-Pr? 

,  P>  Pi 

-  V2P901P1 

irt-. 

Pi  -  Pr 

—  <n 

'Pl-Prr 


P=^{PI-PR)2-r±{PI-Pn 


A  =  |  (P90J  -  Pr, 


)2  +  /3-f  (*W)2. 


(22) 


The  electromechanical  coupling  component  is  given  by 

^es{P,  e)  =  -aeP  -  qeP 2 


where  a  is  the  piezoelectric  coupling  coefficient  and  q  is  the  electrostrictive  coupling  coefficient.  The 
elastic  free-energy  is 

1 

_2 


1, 


Ms)  =  ^Yez 


Figure  7:  1-D  Gibbs  energy  with  three  minima  corresponding  to  the  three  allowed  dipole  states  P_,  P+ 
and  Pgo-  The  1-D  energy  has  four  inflection  points  ±P/,±Pgo/  and  two  local  maxima  (1-D  saddle 
points)  ±PS. 
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where  Y  is  the  Young’s  modulus.  The  total  Helmholtz  free-energy  is  then  given  by 


ip(P,  £)  =  '•PpiP)  +  ^ es{P ,  e)  +  V’ez(e)- 

Balancing  the  internal  energy  ip  and  the  externally  applied  energy  yields  the  Gibbs  free  energy 

G(E,  P,  a,  e)  =  ip(P,  e)  —  EP  —  ae.  (23) 

The  necessary  condition  ^  =  0  yields 

s=Y~1  ( a  +  aP  +  qP 2)  . 

Following  the  same  procedure  as  detailed  in  Section  3,  the  strain  is  substituted  into  the  Gibbs  energy 
(23)  resulting  in 

G(E,  P,  a)  =  $(P,  a)  —  EP  —  ae.  (24) 


4.3.2  Switching  in  the  Absence  of  Thermal  Activation 

Neglecting  thermal  effects,  the  local  polarization  ( P )  for  a  single  crystal  with  an  applied  stress  can  be 
determined  by  solving  the  necessary  conditions 


dG 

dP 


=  0, 


or  by  minimizing  the  Gibbs  energy 


(P)  =  arg  nun  G(E ,  P,  a) 


(25) 


where  the  Gibbs  energy  is  defined  by  equation  (24). 

To  calculate  the  value  of  the  coercive  field,  we  apply  the  necessary  condition  §p  =  0  and  solve  for 
the  value  of  the  electric  field  that  results  in  P  =  Pj.  This  yields 


Ec(a) 


2qPI-af  ,  ,  2 q2Pf 
-^{a)  +  — 


3  aqPf  a2P / 
Y  +  Y 


+  V  (Pr  -  Pi) . 


(26) 


It  is  clear  from  (26)  that  the  coercive  field  is  stress  dependent.  The  coercive  stress  ac  can  be  computed 
by  setting  Ec  =  0  in  (26).  This  results  in 


C  =  - _  \qPi  [nY(Pn  -  Pr)  +  2q2{Pjf  -  ?>aq{PI)2  +  a2(Pr )]  . 


(27) 


The  electric  field  required  to  eliminate  the  90°  minima  in  the  Gibbs  energy  is  denoted  by  E®°.  An 
explicit  expression  for  the  value  of  E^°  is  found  to  be 


El\a) 


-(a) 


a  +  2qPgo 
Y 


2q2Pjo  +  3 agPoo  +  ( a 2  -  mY)p9 o 
Y 


(28) 


As  in  the  2-D  model,  a  direct  180°  ferroelectric  polarization  switch  (i.e.,  P±  — >  P^)  can  occur  if 
E^°  <  Ec.  A  90°  ferroelectric  polarization  switch  (i.e.,  P±  — >  P90)  can  occur  if  E®°  >  Ec  and  a  90° 
ferroelastic  polarization  switch  will  occur  when  a  >  ac. 
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4.3.3  Switching  in  the  Presence  of  Thermal  Activation 


In  the  manner  detailed  in  Section  4.1.3,  thermal  activation  processes  are  incorporated  through  the 
use  of  the  Boltzmann  relation  (14).  For  the  1-D  model,  the  polarization  has  three  allowed  dipole 
states  P-,  P+  and  P90.  The  fraction  of  dipoles  in  each  allowed  orientation  is  given  by  x-,xgg  and  x+. 
Conservation  of  the  number  of  dipoles  yields 


am  +  xgg  +  x+  =  1. 

The  local  polarization  (P)  is  given  by 

(P)  =  X—  ( P- )  +  x90  (P90)  +  x+  (P+)  (29) 

where  (P~) ,  (P90)  and  (P+)  are  the  expected  polarization  values  associated  with  each  allowed  dipole 
orientation.  As  illustrated  in  [24,  31],  these  values  are  found  by  integrating  the  product  of  the  polar¬ 
ization  P  and  the  Boltzmann  probability  density  p(G(P ))  over  the  allowed  polarization  states.  This 
simplifies  to  the  relations 

rPl  Pe~^GdP  (P®0/  Pe~lGdP  fp  Pe~^GdP 

/  p  \  _  J-oo _  /  p  \  _  j--T9o i _  /  p  \  _  J y i 

n-Pr  ’  (-*90/  „ponr  5  (-M-/ 


j~ne~^GdP  '  '  fP9p01  e-^GdP  ’  fp  e-^GdP 


where  1  =  pp-  For  computational  efficiency,  we  can  approximate  these  values  by  use  of  the  necessary 
condition  fp  =  0  which  yields 

(P~)  = - Pr,  (P90)  =  (P+)  = - b  Pr- 

V  V2  r) 

In  this  case,  (P-) ,  (P90)  and  (P+)  are  the  location  of  the  minima  of  the  Gibbs  energy. 

The  time  evolution  of  the  dipole  phase  fractions  are  governed  by  the  first  order  ODE  system 


-P-  P90-  0 

P-  -  (P90-  +P90+)  P+ 

0  p90+  ~p+ 

The  system  results  from  the  assumption  that  transitions  between  the  three  allowed  states  occur  only 
to  the  nearest  neighbor.  The  likelihood  to  switch  out  of  the  P_  orientation  into  the  P90  orientation  is 
denoted  by  p-  and  the  notation  for  the  remaining  likelihoods  is  summarized  in  Table  2. 

The  likelihood  p-  is  calculated  by 

e_ 7  Gi-Pj) 


1 


P-  = 


T  f-D  e-^G(P)dp 

J  —  OO 


where  the  relaxation  time  r  is  the  reciprocal  of  the  frequency  at  which  dipoles  attempt  a  switch.  The 
likelihood  of  switching  out  of  the  P90  orientation  into  the  P_  orientation  is  specified  by  pgg~  and  is 
calculated  by 

\  g-7G(-Pgo) 
m~  ~  x  j^e-iG(P)dp- 

The  likelihoods  pg 0+  and  are  obtained  in  a  similar  manner.  The  likelihoods  can  also  be  evaluated 
in  terms  of  error  functions  as  detailed  in  Appendix  B  resulting  in, 


1 


d-tG(Pj) 


1 


0-a{PI+b/  2)2 


r  /“  e-7 G(p)dP  t  •  erfc  (yfa  (Pi  +  6/2)) 


P+  = 


'yji  (Pr  —  P) 

where  a  =  — -  and  6  =  —2- — - .  The  likelihoods  pgo+,pgo-  and  p-  follow  similarly. 

2  r/ 
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State  Transition 

notation 

P_  - 

4  -P90 

P- 

P90  - 

->  IE 

P90— 

-P90  - 

-  I\ 

P90+ 

P+  - 

4  P90 

P+ 

Table  2:  1-D  transition  likelihoods. 

5  Macroscopic  Polarization  Model 

Nonuniformities  in  the  lattice  structure  due  to  polycrystallinity,  material  nonhomogeneities  and  vari¬ 
ations  across  grain  boundaries  produce  a  distribution  of  Helmholtz  and  Gibbs  energy  profiles  which 
can  be  manifested  as  variations  in  the  local  coercive  held  and  local  remanent  polarization.  Other 
variations  can  be  produced  by  stress  nonhomogeneities  and  variable  effective  fields. 

To  incorporate  these  effects  on  route  to  constructing  a  macroscopic  model,  we  consider  the  coercive 
held  Ec  to  be  a  manifestation  of  an  underlying  density  v\  (Ec)  rather  than  hxed  values  which  is  typically 
assumed  for  single  crystals  having  a  uniform  lattice  structure.  To  create  a  macroscopic  model  for  the 
polarization,  we  also  consider  the  variation  of  effective  helds  in  the  material.  As  detailed  in  [24,  31], 
an  applied  held  E  in  a  ferroelectric  material  is  augmented  by  an  interaction  held  Ej  generated  by 
neighboring  dipoles  which  produce  nonhomogeneous  effective  helds  in  the  material.  This,  along  with 
various  other  processes,  produces  variations  in  the  applied  held  that  can  signihcantly  alter  the  resulting 
polarization.  To  incorporate  these  variations,  we  consider  the  effective  held  Ee  =  E  +  Ej  to  be 
distributed  about  the  applied  held  E  with  an  underlying  density  for  the  interaction  held  Ej  which  we 
denote  by  z^( Ej ).  The  introduction  of  variations  in  the  effective  held  produces  domain  switching  in 
advance  of  the  remanence  point  in  accordance  with  observations  from  experimental  data. 

The  complete  macroscopic  polarization  model  for  nonhomogeneous,  polycrystalline  materials  with 
distributed  coercive  and  effective  helds  is 

/*00  /‘OO 

P(E)  =  /  /  <P  ){E-EI,Ec)u1{Ec)iy2{EI)dEIdEc  (30) 

J  0  J — oo 

where  (P)  is  the  local  polarization  kernel  given  by  (6), (12), (19)  or  (25)  when  thermal  effects  are 
negligible  or  by  (16)  or  (29)  when  incorporating  thermal  relaxation. 


6  Model  Validation  and  Properties 


In  this  section,  the  ferroelastic  models  presented  in  Sections  4.2  and  4.3  are  compared  with  experi¬ 
mental  PLZT  data  reported  in  [13].  Both  ferroelastic  switching  models  employ  the  thermal  activation 
models  presented  in  Sections  4.1.3  and  4.3.3  to  evolve  the  polarization  in  time.  However,  since  the 
measured  values  exhibit  minimal  relaxation,  both  the  1-D  and  2-D  switching  models  were  computed 
with  a  large  7  and  a  small  r  which  resulted  in  negligable  thermal  activation.  As  reported  in  [13] ,  the 
PLZT  composition  is  also  very  near  the  rhombohedral-tetragonal  morphotropic  boundary.  To  model 
this  nonhomogenous  composition,  we  utilize  the  macroscopic  polarization  model  of  Section  5. 

For  the  macroscopic  models,  we  make  the  a  priori  assumption  that  the  density  for  the  coercieve 
held  is  given  by 


vi{Ec(a)) 


C\e ~Ec(a)/bl  e-E1{<r)lbl  >  Q 

0  e-El{<r)/bi  <  o 


(31) 
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and  the  density  for  the  interaction  field  is  given  by 

v2(EI)  =  c2e-El/b*-  (32) 

both  are  plotted  in  Figure  8.  The  relation  Ec(a)  is  specified  by  the  single  crystal  values  (20)  and  (26) 
and  the  variance  of  the  coercive  field  was  estimated  from  the  variance  in  the  experimentally  measured 
values  in  Table  3.  The  variance  of  the  effective  field  was  estimated  by  the  degree  to  which  switching 
occurs  before  remanence  in  the  PLZT  data.  The  parameters  of  the  densities  v\ (Ec)  and  u2(Ei)  used 
in  the  model  were  c±  =  27.0,  b\  =  6.5  x  10_4(MV/m)2,  c2  =  6.0  and  b2  =  5.0  x  10_3(MV/m)2.  To 
implement  the  macroscopic  polarization  given  by  (30),  a  composite  Gaussian  quadrature  was  employed. 

We  note  that  additional  accuracy  can  be  obtained  if  general  densities  u\  and  u2  are  identified  in  the 
manner  detailed  in  [27]  for  180°  field-induced  switching.  Extension  of  the  framework  in  this  direction 
is  under  present  investigation. 

The  ferroelastic  models  were  programmed  in  MATLAB  and  run  on  a  1.7  Ghz  laptop  with  512  MB 
of  RAM.  The  average  computational  times  for  each  model  are  given  in  Table  5.  The  2-D  ferroelastic 
model  requires  the  values  of  the  four  minima  and  the  saddle  points  of  the  Gibbs  energy.  These  nonlinear 
problems  were  solved  using  Newton’s  method  with  a  line  search.  The  computational  cost  of  the  2-D 
model  includes  eight  nonlinear  solves  per  time  step  and  1  LU  factorization  for  the  backward  Euler 
time  discretization  of  the  thermal  activation  models.  The  1-D  ferroelastic  single-crystal  model  does 
not  require  the  nonlinear  Newton  solutions  which  significantly  decreases  the  computational  effort.  It 
is  also  noted  that  when  operating  in  the  absence  of  thermal  activation,  the  local  minimization  models 
given  by  (19)  and  (25)  are  more  computationally  efficient  than  the  thermal  evolution  models  and 
should  be  used  when  appropriate. 

6.1  2-D  Piecewise  Ferroelastic  Model 

The  ferroelectric  hysteron  or  kernel  (12)  is  the  starting  point  for  choosing  several  of  the  ferroelastic 
model  parameters.  The  slope  of  the  hysteron  in  the  180°  switching  regime  is  given  by  r/^1  and  the 
slope  of  the  hysteron  in  the  90°  regime  is  given  by  r/j- 1 .  The  remanent  polarizations  in  the  1  and  3 
direction  are  defined  by  the  parameters  Pin  and  Pm  which  may  have  different  values  depending  upon 
the  internal  structure  of  the  material.  The  inflection  point  P31  is  chosen  to  determine  the  coercive  field 
via  relation  (13).  The  remaining  coefficients  <23,(73  and  I3  affect  the  slope  and  intercept  of  the  linear 


Figure  8:  Typical  Macroscopic  polarization  model  densities  v2(Ei)  and  vi(Ec).  (a)  Normally  dis¬ 
tributed  interaction  field  v2(Ej)  and  (b)  normally  distributed  coercive  field  v\{Ec)  for  varying  applied 
stresses.  Here  Ec( 0)  =  0.35,  Ec(— 6)  =  0.27,  Ec(— 10)  =  0.22,  and  Ec{— 15)  =  0.15. 
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<73  (MPa) 

£m  (C/m2) 

Ec  (MV/m) 

£9°  (MV/m) 

0 

0.247 

0.35  ±  0.01 

0.35  ±  0.01 

-6 

0.247 

0.26  ±  0.01 

0.35  ±  0.01 

-10 

0.235 

0.19  ±  0.02 

0.35  ±  0.02 

-15 

0.215 

0.15  ±  0.02 

0.39  ±  0.02 

Table  3:  Experimental  values  for  03,  P3r,  E®° ,  and  Ec  for  rhombohedral  PLZT  from  data  reported 
in  [13]. 


Ec(a)  relation  (20).  These  values  can  be  ascertained  by  fitting  the  Ec(a)  relation  to  experimental 
data  as  illustrated  in  Figure  9.  Finally,  P\j  is  chosen  so  that  the  90°  switching  occurs  appropriately 
as  the  stress  is  applied.  This  can  be  achieved  by  plotting  the  linear  E^°(a)  relation  to  ensure  that  its 
value  is  sufficiently  large  to  induce  switching  at  the  appropriate  stress  levels.  The  values  of  (/13,  q^i,  a\s 
and  031  can  be  chosen  from  their  effect  on  the  £3  —  £3  relations.  The  ferroelastic  parameters  used  to 
characterize  the  PLZT  data  using  the  2-D  model  are  summarized  in  Table  4. 

The  2-D  model  is  run  with  varying  applied  stresses  in  the  3  direction  as  well  as  an  oscillating 
electric  field  also  in  the  3  direction.  While  the  2-D  model  is  able  to  employ  stresses  and  fields  in 
the  1  direction,  they  are  set  to  zero  to  match  experimental  conditions.  The  reader  is  referred  to  [13] 
for  details  regarding  the  experimental  procedures.  The  behavior  predicted  by  the  2-D  model  (30), 
employing  the  density  choices  (31)  and  (32)  and  with  negligible  thermal  relaxation  is  compared  to 
the  PLZT  data  in  Figures  10  -  12.  It  is  noted  that  a  possible  source  of  error  in  the  model  fit  may 
arise  from  the  rhombohedral  nature  of  the  PLZT  data  since  the  ferroelastic  model  is  derived  in  the 
tetragonal  phase.  The  model  parameters  were  chosen  to  optimize  the  fit  of  the  £3  —  £3  and  £3  —  £3 
data  shown  in  Figures  10  and  11.  The  ferroelastic  model  characterizes  the  90°  switching  that  occurs 
in  the  £3  —  £3  data  as  a  compressive  stress  is  applied  to  the  PLZT  sample.  As  noted  in  Figure  11, 
the  model  also  predicts  a  negative  strain  due  to  an  applied  compressive  stress  as  well  as  the  butterfly 
nature  of  the  £3  —  £3  data.  The  datafit  of  the  <73  —  £3  relation  shown  in  Figure  12  can  be  optimized 
by  setting  the  Young’s  modulus  parameters  to  match  the  slope  of  the  appropriate  part  of  the  173  —  £3 
curve.  However,  this  results  in  an  underprediction  of  the  strains  in  the  £3  —  £3  curve.  Simultaneously 
optimizing  both  the  173  —  £3  and  £3  —  £3  relations  may  be  accomplished  by  using  higher-order  terms 
in  the  electromechanical  coupling  energy  given  by  (1). 

The  PLZT  data  in  [13]  illustrates  the  stress-dependence  of  the  remanence  and  saturation  polariza¬ 
tion.  For  applied  compressive  stresses  greater  than  15  MPa,  the  remanence  and  saturation  polarization 


parameter 

value 

parameter 

value 

P\R 

0.22  (C/m2) 

a3 

0.01  (MV/rn) 

P11 

0.20  (C/m2) 

031 

0.01  (MV/rn) 

P3R 

0.24  (C/m2) 

Ol 

0.01  (MV/rn) 

P31 

0.239  (C/m2) 

013 

0.01  (MV/rn) 

m 

20.0 

13 

500.0  (MV/rn)2 

m 

20.0 

131 

-100.0  (MV/rn)2 

c 

20.0 

qi 

500.0  (MV/rn)2 

Vi 

50000  (MPa) 

113 

-100.0  (MV/rn)2 

>3 

20000  (MPa) 

Table  4:  Parameters  for  2-D  piecewise  ferroelastic  switching  model. 
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a  vs  E 


Figure  9:  03  vs  Ec  for  the  1-D  (solid  line)  and  2-D  (dotted  line)  ferroelastic  models  where  Ec  is  defined 
by  equations  (20)  or  (26). 

significantly  decrease.  Since  the  remanence  polarization  is  a  fixed  model  parameter,  the  present  for¬ 
mulation  of  the  ferroelastic  model  should  be  limited  to  moderate  stress  regimes. 

6.2  1-D  Piecewise  Ferroelastic  Model 

The  stress-free  hysteron  or  kernel  (25)  from  the  ferroelastic  model  provides  the  starting  point  for 
choosing  several  of  the  1-D  ferroelastic  model  parameters.  The  slope  of  the  hysteron  in  the  180° 
regime  is  given  by  ip1  and  the  slope  of  the  hysteron  in  the  90°  regime  is  given  by  rj^ 1 .  The  remanent 
polarization  is  defined  by  the  parameter  Pr  and  the  inflection  point  Pi  is  chosen  to  determine  the 
coercive  field  by  equation  (13).  The  remaining  coefficients  a,  q  and  Y  affect  the  slope  and  intercept 
of  the  linear  Ec(a)  relation  (26).  These  values  can  be  ascertained  by  fitting  the  Ec{cr)  relation  to 
experimental  data  as  illustrated  in  Figure  9.  Finally,  FW  is  chosen  so  that  the  90°  switching  occurs 
appropriately  as  the  stress  is  applied.  This  can  be  done  by  plotting  the  linear  E®°(a)  relation  to 
ensure  that  its  value  is  large  enough  to  begin  switching  at  the  appropriate  stress  levels.  The  values 
of  Pm,  r/i ,  fj  and  A  are  defined  in  equation  (22).  The  parameters  used  to  model  the  PLZT  data  are 
compiled  in  Table  6. 

The  1-D  ferroelastic  model  is  compared  to  the  PLZT  data  in  Figures  13  -  15.  The  model  parameters 
were  chosen  to  optimize  the  fit  of  the  E  —  P  and  E  —  e  data  shown  in  Figures  13  and  14.  As  with  the 
2-D  model,  the  1-D  ferroelastic  model  characterizes  the  90°  switching  that  occurs  in  the  E  —  P  data 
as  a  compressive  stress  is  applied  to  the  PLZT  sample,  but  with  less  computational  cost  than  the  2-D 
model.  The  model  also  predicts  a  negative  strain  due  to  an  applied  compressive  stress  as  well  as  the 
butterfly  nature  of  the  E  —  e  data. 

We  note  that  for  this  characterization  regime,  the  simplified  1-D  model  provides  approximately 
the  same  accuracy  as  the  2-D  model  but  is  nearly  2  orders  of  magnitude  faster  as  indicated  in  Table  5. 


Single  Crystal 

Macroscopic 

1-D 

0.3  sec 

67.8  sec 

2-D 

4.6  sec 

1155.4  sec 

Table  5:  Average  computational  times  for  the  ferroelastic  switching  models  running  on  a  1.7Ghz  laptop 
with  512MB  of  RAM.  The  applied  electric  field  was  E  =  £0  sin  (37rf )  where  t=[0,l]  with  At  =  0.005. 
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parameter 

value 

parameter 

value 

Pr 

0.24  (C/m2) 

a 

0.01  (MV/rn) 

Pi 

0.239  (C/m2) 

q 

500.0  (MV/rn)2 

P901 

0.015  (C/m2) 

Y 

20000  (MPa) 

V 

20.0 

m 

20.0 

Table  6:  Parameters  for  the  1-D  piecewise  ferroelastic  switching  model. 


g3  =  0  MPa  o3  =  -6  MPa 


(a) 


(b) 


-0.6  -0.4  -0.2  0  0.2  0.4  0.6 

Eg  (MV/m) 


-0.6  -0.4  -0.2  0  0.2  0.4  0.6 

E3  (MV/m) 


(c)  (d) 

Figure  10:  Macroscopic  £3  —  £3  relation  for  varying  applied  stresses.  Homogenized  2-D  ferroelastic 

model  ( - )  and  experimental  PLZT  data  (•  •  • ):  (a)  cr3  =  0  MPa,  (b)  cr 3  =  —6  MPa,  (c)  <73  =  —10 

MPa,  and  (d)  <73  =  —15  MPa. 
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E3  (MV/m)  E3  (MV/m) 

(a)  (b) 

c3  =  -10  MPa  o3  =  -15  MPa 


_0.6  -0.4  -0.2  0  0.2  0.4  0.6  -0.6  -0.4  -0.2  0  0.2  0.4  0.6 

E3  (MV/m)  E3  (MV/m) 

(c)  (d) 


Figure  11:  Macroscopic  E3  —  £3  relation  for  varying  applied  stresses.  Homogenized  2-D  ferroelastic 

model  ( - )  and  experimental  PLZT  data  (•  •  • ):  (a)  <73  =  0  MPa,  (b)  <73  =  —6  MPa,  (c)  <73  =  —10 

MPa,  and  (d)  <7 3  =  —15  MPa. 


Figure  12:  (a)  Macroscopic  <73  —  P3  relation  and  (b)  macroscopic  03  —  £3  relation.  Homogenized  2-D 
ferroelastic  model  ( - )  and  experimental  PLZT  data  (•••). 


Figure  13:  Macroscopic  E  —  P  relation  for  varying  applied  stresses.  Homogenized  1-D  ferroelastic 

model  ( - )  and  experimental  PLZT  data  (•  •  • ):  (a)  <73  =  0  MPa,  (b)  <73  =  —6  MPa,  (c)  <73  =  —10 

MPa,  and  (d)  03  =  —15  MPa. 
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-0.6  -0.4  -0.2  0  0.2  0.4  0.6 

E3  (MV/m) 


-0.6  -0.4  -0.2  0  0.2  0.4  0.6 

E3  (MV/m) 


(c) 


(d) 


Figure  14:  Macroscopic  E  —  £  relation  for  varying  applied  stresses.  Homogenized  1-D  ferroelastic 

model  ( - )  and  experimental  PLZT  data  (•  •  • ):  (a)  <73  =  0  MPa,  (b)  <73  =  —6  MPa,  (c)  <73  =  —10 

MPa,  and  (d)  <73  =  —15  MPa. 


°3vsP3  a3vse3 


(a)  (b) 

Figure  15:  (a)  Macroscopic  a  —  P  relation  and  (b)  macroscopic  a  —  e  relation.  Homogenized  1-D 
ferroelastic  model  ( - )  and  experimental  PLZT  data  (•••). 
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7  Concluding  Remarks 


This  paper  addresses  the  development  of  a  framework  for  characterizing  stress-induced  90°  and  180° 
switching  inherent  to  ferroelectric  materials  in  a  manner  suitable  for  subsequent  transducer  and  control 
design.  The  model  builds  upon  the  homogenized  energy  framework  developed  for  ferroelectric  materials 
in  [31]  and  extended  to  provide  a  unified  characterization  framework  for  ferroic  compounds  in  [29,  30]. 

In  the  first  step  of  the  development,  three  techniques  are  employed  to  construct  Helmholtz  and 
Gibbs  energy  functionals  used  to  characterize  the  electromechanical  behavior  of  homogeneous,  single 
crystal  compounds.  The  first  exploits  a  classical  Landau-Devonshire  energy  formulation  whereas  the 
second  employs  a  piecewise  quadratic  formulation  to  eliminate  numerical  and  physical  restrictions 
inherent  to  high-order  polynomials  comprising  the  Landau  Devonshire  functional.  Both  choices  quan¬ 
tify  the  2-D  polarization  P  as  a  function  of  2-D  stress  and  field  inputs.  The  third  choice  is  similar  to 
that  employed  for  SMA  [14,  23,  24,  29]  in  the  sense  that  it  is  1-D  with  three  wells  corresponding  to 
±180°  and  90°  equilibria.  The  construction  of  this  functional  is  phenomenological  but  the  resulting 
decrease  in  dimension  significantly  diminishes  implementation  time.  In  all  three  cases,  the  functionals 
are  directly  minimized  to  provide  kernels  for  characterization  in  the  absence  of  thermal  relaxation 
or  balanced  with  the  relative  thermal  energy  through  Boltzmann  principles  to  incorporate  relaxation 
phenomena. 

In  the  second  step  of  the  development,  the  effects  of  material  nonhomogeneities,  polycrystallinity, 
and  variable  effective  fields  are  incorporated  by  assuming  that  properties  such  as  local  coercive  and 
interaction  fields  are  manifestations  of  underlying  distributions  rather  than  constants.  Stochastic 
homogenization  in  this  manner  provides  a  macroscopic  model  which  characterizes  field  and  stress- 
induced  90°  and  180°  switching  in  ferroelectric  materials. 

In  the  present  formulation,  a  priori  choices  of  normal  densities  are  made  when  quantifying  lo¬ 
cal  coercive  and  interaction  field  distributions.  This  yields  macroscopic  models  with  a  small  number 
of  parameters  to  be  identified  but  restricts  the  accuracy  of  the  framework  to  realizations  of  the  a 
priori  choices.  Alternatively,  one  can  employ  general  densities  to  be  identified  through  least  squares 
techniques  in  the  manner  detailed  for  field-induced  180°  ferroelectric  switching  in  [27].  This  greatly  in¬ 
creases  the  flexibility  and  accuracy  of  the  model  while  maintaining  the  same  complexity  with  regard  to 
subsequent  control  design.  Extension  of  the  framework  in  this  direction  is  under  present  investigation. 

One  way  in  which  the  homogenized  energy  framework  has  been  employed  for  control  design  is 
through  the  construction  of  inverse  filters  which  can  be  used  to  linearize  the  transducer  response. 
Details  illustrating  this  technique  in  the  context  of  the  magnetic  model  quantifying  180°  field-induced 
switching  can  be  found  in  [20]  and  development  of  inverse  filters  for  the  present  theory  constitute  a 
second  facet  of  current  research. 


Appendix  A:  Coefficients  of  Stress  Dependent  Gibbs  Energy 

We  summare  here  the  coefficients  employed  in  the  energy  relations  of  Section  4.2.1.  The  coefficients 
result  from  substituting  the  equations  for  the  strains  into  the  Gibbs  energy  relation  and  results  in  the 
following: 
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Appendix  B:  Conversion  of  Polarization  Integrals  to  Error  Functions 


For  fixed  temperatures,  we  employ  the  piecewise  quadradic  relation  Helmholtz  energy 


V>(P) 


'  Iv(P  +  Pr)2 
<  §*/(P-P*)2 

k  \ti(Pi-Pr)(^-Pr) 


P  <  -Pj 
P>  Pi 
\P\  <  Pi 


(34) 


where  Pj  and  Pr  denote  the  positive  inflection  point  and  polarization  at  which  the  minimum  occurs. 
The  resulting  Gibbs  energy  is  derived  by  combining  the  potential  energy  of  a  dipole  in  the  field  with 
the  Helmholtz  energy  throughout  the  lattice  to  yield 


G(P,  E)  =  *1>(P)  -  EP. 


(35) 


In  the  thermal  evolution  model,  integrals  of  the  form 

/“  e-^P’E)dP  and  Per~'G(RE>dP 


(36) 


must  be  computed  where  7  =  \/  Ap-  Since  ^  is  a  piecewise  quadradic,  the  integrals  are  Gaussians 
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and  can  be  represented  as  error  functions  for  computational  efficiency.  It  is  helpful  to  note  that 
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where  a  =  Pp  and  b  =  —2 (Pr  —  p).  Completing  the  square  results  in 
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which  gives 


For  the  next  integral,  it  is  helpful  to  note  that 
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Following  the  procedure  above,  we  arrive  at 
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where  a,f3  and  b  are  previously  defined.  Multiplication  by  unity  followed  by  addition  and  substruction 
of  the  same  expression  yields 
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Thus  the  final  relation  is 
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